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We present a technique for the dissipative preparation of highly entangled multiparticle states of atoms cou¬ 
pled to common oscillator modes. By combining local spontaneous emission with coherent couplings we en¬ 
gineer many-body dissipation that drives the system from an arbitrary initial state into a Greenberger-Horne- 
Zeilinger state. We demonstrate that using our technique, highly entangled steady states can be prepared effi¬ 
ciently in a time that scales polynomially with the system size. Our protocol assumes generic couplings and will 
thus enable the dissipative production of multiparticle entanglement in a wide range of physical systems. As an 
example, we demonstrate the feasibility of our scheme in state-of-the-art trapped-ion systems. 
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Multiparticle entanglement is an essential resource for 
quantum computation and information JT], e.g. in quan¬ 
tum error correction m a, quantum memories El, and 
entanglement-enhanced quantum measurement schemes 0 
a. Among the most important states which exhibit genuine 
multiparticle entanglement are Greenberger-Home-Zeilinger 
(GHZ) states 

|GHZ) = ^ (1000...0)+ 1111...1)). (1) 

Deterministic preparation of such states has so far been per¬ 
formed using time-dependent unitary gates GHM, which 
have recently yielded impressive progress towards entan¬ 
gling larger numbers of qubits ciHia, and feedback con¬ 
trol schemes OlHIbl. These operations however suffer from 
quantum noise, causing decoherence and dissipation so that 
it remains difficult to prepare high-fidelity multiparticle en¬ 
tangled states with these methods Eiiia. Recently, dissi¬ 
pative state preparation has been proposed as an alternative 
approach where the dissipative environment is actively engi¬ 
neered and used to prepare states relevant for quantum infor¬ 
mation and simulation lfT9l - E3l . Numerous theoretical studies 
on the production of bipartite entangled states EHm have 
since been performed and the first experimental demonstra¬ 
tions ll32l - lT5l have been realized. More recently, also dissi¬ 
pative schemes for the generation of multipartite entangled 
states II 36 H 5 TII have become available, e.g. for the preparation 
of states stabilized by local interactions \23\ l42l . It has 
however remained a challenge to prepare in a scalable way 
states, like GHZ, that are highly entangled in the sense that 
they cannot be stabilized by local operators Il2ni23ll40ll . 

In this Letter, we extend the range of the dissipative ap¬ 
proach by demonstrating a scalable technique for the dissi¬ 
pative preparation of highly entangled states of many parti¬ 
cles. We show that, by using local spontaneous emission as 
a generic source of dissipation, we can engineer nontrivial 
many-body dissipative interactions HSl which are tailored to 
produce multiparticle GHZ states. Our scheme is determin¬ 
istic and operates by continuous optical driving from an arbi¬ 
trary initial state towards the desired steady state using weak 
classical fields. The preparation time of our protocol is found 



FIG. 1. Protocol. Preparation of a GHZ state of N qubits is realized 
by two operations: (ij “Z-Pumping” of states with other than 0 or A 
atoms in state |1) to |0)®^. (ii) |0)®^ is a superposition of |GHZ) 
and|GHZ_) = (|0)®^ - |l)®^)/y2 so that |GHZ_) needs to be 
removed by the parity-selective “X-pumping”. 

to exhibit a favorable polynomial scaling with the number of 
qubits. In addition to our generic system-independent scheme, 
we describe an implementation in a system of trapped ions. 

In our protocol, preparation of a steady GHZ state Q of A 
qubits starting from an arbitrary initial state is accomplished 
by two simultaneous operations shown in Fig. [T] f/j Pumping 
all states with more than zero but less than all qubits in state 
|1) (0 < ni < N) to the state |0)®^, which can be writ¬ 
ten as a superposition jO)®'^ = |GHZ) + |GHZ_), and (ii) 
removing the GHZ state with the wrong phase, |GHZ_) = 
(| 0 )®Af _ from the subspace spanned by |0)®'^ 

and |1)®^. Operation (i) is implemented such that it fulhlls 
the main requirement of a dissipative protocol: it has to pump 
an exponential number of states efficiently, i.e. in polynomial 
time. In principle, standard optical pumping ll52l satisfies this 
criterion as well, but would also erase the state |1)®^, thus 
ruling out the possibility to prepare |GHZ) with high fidelity. 
Instead, we design a new procedure which is selective in the 
number of atoms in |1). This allows us to pump only states 
with 0 < ni < N to states with ni — 1, and thus eventually 
to |0)®^ (ni = 0). We refer to this operation shown in Fig. Ih 
as “Z-pumping”, since it is based on counting the number of 
atoms in the eigenbasis of Z = |0)(0| — |1)(1|. 

A second operation (ii) is required to remove the undesired 
|GHZ_) state in a continuous manner, as also illustrated in 
Fig.B To this end, we perform a pumping process selective 
in the parity V = (A, = |1),(0| + |0)a(l|). Here, 
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we apply the recipe for stabilizer pumping from Ref. m 
to the case of the parity stabilizer, V: Expressed in terms of 
the eigenstates of X, |±) = (|0) ± \l))/\/2, any state |'0) 
with 7^|'0) = +1|'(/'), such as |GHZ), is a superposition of 
only those product states that contain an even number, n_, 
of |—) qubits; e.g., for N = 3, |GHZ) = (| + ++) + | + 

-) + I —I—) + 1-1"))/2. On the other hand, n_ is odd 

for any state {ip) with Vlpj) = —1|'0)> such as |GHZ_); for 
N = 3, |GHZ_) = (I + +-) + I + -+) + I - ++) + I - 

-))/2. By pumping all states with odd n_ to other states, 

in the following referred to as “X-pumping”, we achieve the 
depumping of |GHZ_). 

The two operations required for our protocol can be realized 
using a generic setup as described next. We assume a general 
system of N particles (“atoms”), shown in Fig. |^(a)-(c). Each 
atom supports two stable ground states |0) and |1) and two 
excited states |e) and |/). The atoms are driven by classical 
multi-tone driving fields, with identical amplitudes on all ions, 
as described by the Hamiltonians 


ffiF) _ l^(F) 

drive,Z “ 2 ^ ^ 


N 


J2\e)a{M+H.c. ( 2 ) 


a—1 

N 


^ |/),(-| -f H.C., (3) 


a—1 


(p'\ (F) 

with strengths ' and detunings A; where I = Z,X de¬ 
notes the desired operation and Eghz the field tone. The tran¬ 
sitions of the atoms are collectively coupled to two harmonic 
oscillator modes, b and c, e.g., two resolved resonator modes 
in cavity or circuit QED CIlIol, or two phononic modes in an 
ion trap setup EnnHia, as modeled by the Hamiltonians 


N 

Hint,z = 9b'^'^\^)a{e\+H.C., (4) 

a—1 
N 

Hint,X = gc^'^\-)a{f\+ H.C., (5) 

a—1 

with g being the coupling constant. For the dissipative process 
we consider decay by spontaneous emission from the excited 
states j to the ground states i at a rate 7 ^- (with 7 ^ = ^^ 7 ^^), 
described by jump operators L^.-^a = y/^\i)a{j\ for i € 
{ 0 , 1 }, j G {e, /} acting incoherently on each atom a. 

We realize the Z- and X-pumping operations by engineer¬ 
ing selected transitions to be driven resonantly, while sup¬ 
pressing others due to off-resonant driving. For Z-pumping 
we use the coupling configuration in Fig. [^(b). Here, a cou¬ 
pling g of the oscillator b to the transition |e) O | 1 ) and 
a weak drive on the same transition are used to effectively 
“count” the number of atoms in state |1). In Fig. s (d), 
we illustrate this mechanism for N = 3 qubits. The weak 
driving tones ' couple the ground states to atom-excited 
states. For example we consider |'0iio) = |110) which is cou¬ 
pled to IV’iio.e) = (|el 0 ) + |le0))/-\/2- This state is in turn 
coupled to an oscillator-excited state |'0iio)|l)t) by the atom- 
oscillator coupling. Because of constructive interference be¬ 
tween the two terms in |' 0 iio,e) this coupling has a strength 



d) ^(11.> + |1<.1) + |.11)) |^,„,) = J=(l.0) + |d0» 1^00) 



FIG. 2. Setup (a)-(c) and dissipative mechanism (d) for GHZ prepa¬ 
ration, shown for N = 3 qubits. Setup (a). We consider a chain of N 
sub-systems (“atoms”) with four levels, coupled to two common har¬ 
monic oscillators. As dissipative processes we assume spontaneous 
emission from the excited levels 7e//' We apply coupling configura¬ 
tions ‘Z’ (b), and ‘X’ (c) consisting of 2{N — 1) driving tones 
with detunings in (b)> 2 [(A-l-1)/2J tones with ™d aF 
in (c), and couplings of the atoms to the oscillator modes g. (d) “Z- 
pumping” towards |0)®^ using the couplings in (b). Ground states 
are coupled to atom-excited states by weak driving. Depending on 
the number ni of atoms in 11), the excited states form dressed states 
with oscillator-excited states at energies rty/nig. By applying fields 
with detunings = '/Fg for 1 < F < N — 1 (only the red- 

detuned fields are shown), all states except |GHZ) and |GHZ_) are 
pumped to |0)®^ = (|GHZ) -f |GHZ_))/V^. |GHZ_) is emptied 
by the parity-selective X-pumping as described in the text. 


of v/2g. As a consequence of the strong atom-oscillator cou¬ 
pling, the atom- and the oscillator-excited state form dressed 
states IV'iio.i) = (IV’iio^)| 0 )b±|V'iio)|l)h)/^at detunings 
A± = ±-\/2p, see Fig. ^(d). Applying a weak driving field 
with a detuning \Az\ = Fig, one thus excites ground states 
with ni = 2 to excited states like IV'iio.e)- Since |e) is subject 
to spontaneous emission to | 0 ) and | 1 ), |' 0 iio,e) decays either 
back to the manifold of states with ni = 2 or “forward” to 
states with ni = 1. In general, the couplings of the atom- and 
oscillator-excited states have a strength of that depends 
on the number ni of atoms in 11 ) so that the dressed states are 
shifted by A± = ±^/nig. This creates a resonance condition 
depending on ni of the initial state which we can use to selec¬ 
tively drive a manifold. Applying a drive with \Az\ = v/Tg, 
states with ni = 1 are thus pumped to ni = 0 , and thereby 
to I GHZ) and |GHZ_). On the other hand, with the detun¬ 
ings Az = ±5 and -AFlg, the state | 1 )®^ is excited only 
off-resonantly and thus decays slowly (see Fig. [^(d), so that 
I GHZ) remains almost unaffected. 

To realize the full Z-pumping process based on the mecha¬ 
nism above, we apply a weak drive consisting of 2 (A — 1 ) 
tones FF with detunings A^^ = aFf g ranging from 
F = 1 to F = N— 1. This gives rise to effective decay 
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processes ll5^ described by 

4o:i2 = (l<rii<iV-l). (6) 

Here, is the projector onto the ground states with ni 
atoms in state |1), a denotes the atom subject to decay, and 

= 7 oe(H^^^/ 7 e)^ are the strongly enhanced decay 
rates of the states which are resonantly excited. Each ground 
state with 1 < ni < TV — 1 then decays towards one with 
ni — 1 by the engineered spontaneous emission in Eq. 

The concatenation of these decay processes causes a contin¬ 
uous drift towards states with smaller rii, finally ending at 
the state with ni = 0. For a suitable choice of the 

15^ , the preparation time r of |0)®^ and the corre¬ 
sponding rate r+ = l/r can then have a favorable scaling 
T oc l/r_|_ ^ log(iV), which is similar to optical pump¬ 
ing where the transition rate of each level is proportional 
to the number ni of its excitations. However, as opposed 
to standard optical pumping, the Z-pumping is engineered 
such that it does not affect and thus neither |GHZ), 

since resonant excitation out of |1)®^ would require the tone 
F = N. By excluding such tones from the drive the ex¬ 
citation out of |1)®^ is off-resonant and thus much weaker. 
The weak off-resonant excitation of |1)®^ by the other driv¬ 
ing tones leads only to a small leakage from |GHZ). Since 
the energy gap between the resonances and the driving tones 
(■v/TV—yTii )5 decreases with N, the leakage rate from |GHZ) 
increases with the number of qubits, and can be estimated to 
be There- 

suiting error can, however, be compensated by reducing the 
speed of the protocol by a small polynomial factor as we dis¬ 
cuss below. 

Having operation (i) realized, we now turn to operation (ii), 
the X-pumping, which removes states with an odd number n_ 
of atoms in |—). To implement it we use a similar mechanism 
as was used to “count” ni above, except that now we need to 
do this in a different basis. We achieve this operation using 
the coupling configuration in Fig. |^(c): Coupling the transi¬ 
tions from |0) to I/) and from |1) to |/) by fields of the same 
strength, but with the opposite phase, results in a coupling of 
the transition from |—) to |/). To avoid interference with the 
Z-pumping, we consider a second excited state |/) and a sec¬ 
ond oscillator mode c. It is nonetheless possible to implement 
the described operations with a single excited level and a sin¬ 
gle oscillator mode in a stroboscopic manner, resulting in a 
quasi-steady state. 

In the X-pumping, the transitions between the excited level 
I/) and the ground levels |0) and |1) are coherently coupled 
to the harmonic oscillator c and excited by a multi-tone drive 
'. Opposite phases on both transitions result in a coupling 
of the transition |/) o |—), similar to the coupling |e) O |1) 
in the Z configuration. In this way, we make the X-pumping 
selective in n_ in a similar manner as the Z-pumping is selec¬ 
tive in ni: Applying 2 [(TV + 1)/2J field tones with detunings 

= ±\/Fg for F = 1,3,5,... {F < N) we resonantly 
excite |GHZ_) to dressed states which lie at FyPnZg for 
n_ = 1, 3, 5,... (n_ < N). Thereby we make |GHZ_) de¬ 
cay to random states by effective spontaneous emission with 


a strong rate r:|f ~ 2(H^^)^/7/, Similar to the Z-pumping, 
the decreasing energy gap between the dressed states gives 
rise to a leakage rate from GHZ, rf ~ N‘^-^fVl\lg‘^ (using 

= VLx for odd F and = 0 for even F), which 
increases with the number of qubits TV. 

The simultaneous action of the Z- and the X-pumping pre¬ 
pares I GHZ) and maintains it as the unique steady state of the 
dissipative dynamics. However, since the Z-pumping is dis¬ 
turbed by the X-pumping, the latter has to be sufficiently weak 
so that the Z-pumping has a sizable probability of reaching the 
final state |0)®^ before being subject to X-pumping; picking 
similar rates for the X-pumping and the total Z-pumping r+, 
this requirement does not slow down the preparation process 
significantly ll5Tl . To find the preparation time we can con¬ 
sider a simple model where the rate of pumping to |0)®^ 
is determined by T^, and where |GHZ) and |GHZ_) are 
pumped out with rates r_ = + r:5 and E;^, respectively. 

Further details on the model are given in the Supplementary, 
but in short we find that the steady state fidelity F is deter¬ 
mined by the ratio of the decay out due to off-resonant excita¬ 
tion at the rate r_ ~ 'yNil'^/g'^, and the effective preparation 
rate r+ ^ f^^/7 (using ~ fl/y/N log TV, ~ H for 
F < 2TV/3, and ~ ^2{N - F)/Fn for F > 2TV/3). 
This gives a steady state error £ = r_/r_|_ ^ jg^ which 

is approached exponentially in time ^ 

To avoid having an increasing error with the qubit number 
TV we assume that we can control the decay rates of the excited 
states |e) and |/). This is for instance the case if the states are 
metastable states coupled to higher lying unstable states with 
a laser field ll^ l34l . The increase of the error with TV can 
then be compensated by having a sufficiently low decay rate 
7 ~ g\f£ l\fN, which keeps the error £ constant for growing 
TV, but prolongs the necessary preparation time tchz- These 
considerations, however, assume a weak driving H, whereas 
for strong driving the pumping rate becomes limited by power 
broadening. We therefore need to use a suitably low driving 
strength H ^ ^/\fN which sets a limit on the preparation 
rate. These considerations and parameter values can be turned 
into a rigorous upper bound on the preparation time ll5^ : 

tghz oc TV^/^ (logTV) /(sVf). (7) 

We can thus prepare a GHZ state at any desired fidelity 
Tghz = 1 — £ within a preparation time that has a low- 
order (TV^/^) polynomial scaling in the number of qubits, for 
a coupling g independent of TV. If instead the total prepa¬ 
ration time is restricted to t < Tmax, then the prepara¬ 
tion error of an TV-particle GHZ state will necessarily obey 
£ > jg^ and £ > e“Tiiaxr+^ limiting the achievable 
fidelity to Fghz < 1 - Nji^gTmaxf lES. 

To confirm the results of the simple model we simulate the 
protocol numerically; In Fig.[^(a) we plot the time-evolution 
towards a steady GHZ state for TV = 2,..., 8 qubits result¬ 
ing from our protocol. The plots are obtained by numeri¬ 
cally simulating an effective master equation ll54l as well as 
the complete master equation truncated to one or two exci¬ 
tations. Here we optimize the available parameters (driving 
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FIG. 3. Evolution towards a steady GHZ state. Starting from a fully 
mixed state, we numerically solve an effective master equation. The 
curves in (a) show the evolution for two to eight qubits (different col¬ 
ors from blue to black) and are obtained by numerically optimizing 
all available parameters to reach a fidelity of Fghz = 0.9 (black 
dashed line) within as short a time as possible. In (b) we show the 
scaling of the preparation time with the number of qubits. Both pan¬ 
els (a) and (b) show different degrees of truncation of the Hilbert 
space (dash-dots/blue squares - effective ground state dynamics af¬ 
ter adiabatic elimination, solid lines/green circles - one excitation, 
dashes/red triangles - two excitations). In (b), small symbols rep¬ 
resent analytically optimized parameters and large symbols numer¬ 
ically optimized parameters. We find a polynomial scaling of the 
preparation time which is within our analytical bound (black dashed 
line in (b)). 


strengths, tunable decay rates) to reach a fidelity Fqhz = 0.9 
of the GHZ state in minimum time. The resulting prepara¬ 
tion times and the analytical bound from Eq. (0 is shown in 
Fig. [^(b). These results confirm that the scheme exhibits a 
low-order polynomial scaling of the preparation time with the 
number of qubits. In contrast to the scheme in ll40l . our proto¬ 
col requires only two operations for a GHZ state of N qubits. 
Furthermore, the highly directed Z-pumping is only weakly 
perturbed by the polynomially weaker X-pumping and thus 
allows for a polynomial scaling of the protocol. 

The ingredients necessary for our scheme are available in 
trapped ion experiments. One suitable setup consists of a 
chain of N trapped ions, each with two (meta-)stable ground 
levels |0) and |1) and two auxiliary levels, |e) and |/). Tun¬ 
able decay of the auxiliary levels by spontaneous emission 
can be realized by a repumper to a higher lying rapidly de¬ 
caying state Ealii. Two phononic modes, cooled to the 
ground state, and coupled to the sidebands of the ions’ tran¬ 
sitions, can be used as the harmonic oscillators b and c. For 


the pumping, we require 2{N — 1) field tones in the Z con¬ 
figuration and 2[(A^ + 1)/2J tones in the X configuration. 
An alternative stroboscopic implementation requires only a 
single auxiliary level, interchanging between the roles of |e) 
and I/), a single phononic mode, interchanging between be¬ 
ing b and c, and a single field tone with tunable detuning ap¬ 
plied on either the transition |1) o |e) or |—) -o- |/). With 
5/27r ~ 10 kHz, typical preparation times are r ~ 30 ms. On 
such timescales collective dephasing needs to be taken into 
account, but can be overcome by switching the role of |0) and 
|1) in the Z-pumping for every second ion, thereby preparing 
(|0101...) -f |1010...))/-\/2 which is in a decoherence free sub¬ 
space and is equivalent to a GHZ state ini, or by using clock 
states m- Fluctuations of g of 1% result in a reduction of 
the fidelity by 0.01 — 0.1 for N = 2,8 qubits, whereas 
fluctuations of 0.1% have an effect at the sub-percent level. 
AC Stark shift fluctuations are suppressed since both red- and 
blue-detuned driving tones (e.g. A± = ±yTfi5) are applied. 
Heating of the motion would constitute another error on the 
timescale of the scheme, but this can be made negligible for 
cryogenic traps Ea. 

We have shown that dissipative state preparation can be ex¬ 
tended to the efficient generation of highly entangled steady 
states of many particles. We achieve this by engineering com¬ 
plex multiparticle dissipation which deterministically drives 
the system into a desired steady state within a time scaling 
only polynomially with the size of the system. The generic 
couplings assumed in our approach can be found in a vari¬ 
ety of physical systems, such as trapped ions where the ba¬ 
sic ingredients of the scheme have already been demonstrated 
ll34ll . As specific examples we have considered the prepara¬ 
tion of highly entangled GHZ states, which are paradigmatic 
multiparticle entangled states, but the developed techniques 
are applicable to a range of other quantum information tasks. 
Particularly relevant further possibilities are the construction 
of quantum etTor con'ecting codes Elia or the observation of 
exotic phase transitions ll22l induced by multiparticle dissipa¬ 
tion. 
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SUPPLEMENTAL MATERIAL “SCALABLE DISSIPATIVE PREPARATION OP MANY-BODY ENTANGLEMENT” 

In this Supplemental Material to the Letter “Scalable dissipative preparation of many-body entanglement” we present the dy¬ 
namical model of the generic system under consideration (Section [I|, the coupling conhgurations for the described protocols 
(Sectionjl^ and the resulting effective dynamics of the system (Sectionjl^. We discuss the engineering of the dissipative many- 
body interactions for GHZ state preparation (Section [TV]) and strong driving effects (Section |V]i. In Section|V^we analyze the 
scaling of the protocol both for weak and for strong driving. 
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I. DYNAMICAL MODEL OF THE SYSTEM 
The dynamics of the open system is modeled by a master equation of Lindblad form 

P = ^i.P) = -* [H, P] + X! ^kP^k - \ {^iLkP + pLlLk^ . (S8) 

k 

Here, p is the density matrix of the system and C denotes a Liouvillian of Lindblad form. The Hamiltonian H contains the 
coherent interactions, while sources of dissipation present in the system are represented as Lindblad (‘jump’) operators L^- The 
Hamiltonian for the coupling configurations detailed in Fig. 2 b)-c) is, in its most general form, given by 

H = + iTint + ^^drive- (S9) 

We consider a free Hamiltonian i7free which contains the energies of the levels of the N atoms and the harmonic oscillator 
modes, 

iFfree = <^eJee + f Jf f + + WcC^C. (SIO) 

Here we have introduced Jij = I*)a0l and made the simplifications wq = wi = 0 and h = 1. An 

interaction Hamiltonian describes the atom-oscillator coupling, and a drive Hamiltonian TTdrive contains the fields used to per¬ 
form coherent excitations of the system. The interaction terms and drives required for GHZ preparation are detailed in Section]^ 

The excited degrees of freedom in the system are generally subject to dissipation. Here the excited states |e) and |/) of each 
atom undergo spontaneous emission to each of the ground states |0) and |1), described by the jump operators 

L^o.,a = y/l^e\^)a{e\ (Sll) 

^Tie.a = V^|l)a(e| (S12) 

^70/.a = VW|0)a(/| (S13) 

= V^|l)a(/|, (S14) 

where the subscript a denotes the atom number. The total decay rates of the excited levels are given by 7 e = 7oe + 7ie and 
7 / = 7o/ + 7i/- For simplicity we will assume equal decay rates to both ground states, which is, however, not crucial for the 
protocol. In addition, we assume decay of excitations of the two oscillator modes, b and c, represented by 

(S15) 

Lk, = y/i^cC. (S16) 

The scheme does not require oscillator decay at all {Kb = Kc = 0). It may, however, still be useful to avoid heating. 


II. COUPLING CONFIGURATIONS 


The atom-oscillator couplings of the four coupling configurations in Fig. 2 b)-c) of the main text are described by the interaction 
Hamiltonian 


-^int — Hint,Z H" ^int^X 

(S17) 

~ 9 Jle “1“ 

(S18) 

Hint,x = g + bJif'j 

(S19) 

Here, ‘Z’ and ‘X’ denote the coupling configurations introduced in the main text (recall that |—) = (|0) — |l))/v^- These terms 
describe that an atomic excitation can be exchanged coherently with the respective harmonic oscillator with a coupling constant 
g. By iTint.z an atomic excitation in |e) is exchanged with the oscillator b, leaving the atom in |1). iFint.x couples the excited 
level 1/) to |—) = ;( 72 (I^) + 1^)) while exchanging the atomic excitation with the oscillator c. The coherent excitation of the 
atoms by classical driving fields is modeled by the drive Hamiltonian 

-^drive ^drive,Z -^drive,X 

(S20) 

iFdrive.Z = ^ ^ Jel + H.C. 

F 

(S21) 

iFdrive.X = ^ ^ + H.C. 

F 

(S22) 
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(ip\ (F) 

Here, we generally allow for several field tones with Rabi frequencies ' and frequencies from which we will deduce 

( F) 

detunings ' of the respective fields. 


III. EFFECTIVE DYNAMICS OF THE OPEN SYSTEM 

In the following, we provide a detailed analysis of the effective dynamics of the system described in Section To this end, 
we briefly introduce the effective formalism presented in Ref ll5^ and use it to derive effective operators for the coupling 
configurations presented in Section [11] 


A. Effective operator formalism 


As can be seen from Section [I] the dissipation affects the excited levels |e) and |/) and the oscillator modes a and b. For 
weak driving the decaying degrees of freedom can be adiabatically eliminated from the master equation. This is done using the 
effective operator formalism presented in Ref. Il54l . In this way, the dynamics of the master equation are reduced to effective 
couplings between the ground states of the system, described by an effective master equation 


P — '^eff(p) — t [Tfeff 5 p] T ^ ^ 2 ( “f ■ (S23) 

k 

Since we are dealing with multiple field tones F that give rise to the effective couplings, we use the extended formalism for 
many fields (cf. El) with 


= + H.c. (S24) 

F 

ife.eff (S25) 

F 

= iTfree + iFint - 


/ C'\ f _ f C'\ 

Here, ' denotes the exciting part and V_ ^ the deexciting part of the drive V± = ■ The non-Hermitian Hamiltonian 

which contains the frequency u!^^\ describes the time evolution of the excited states. The drives are regarded as perturba¬ 
tions and thus denoted by “I/”; they are defined by the drive Hamiltonians iTdrive in Sectionj^so that we use V = Tfdrive when 
we derive the effective operators for each of the coupling configurations below. Based on the assumption of weak excitation we 
will restrict our discussion to the states which have at most one atomic or oscillator excitation. I/nh then contains the energies 
and couplings of the excited states. Since needs to be inverted to compute the effective operators, we will, for each coupling 
situation, start out by discussing this entity. 


B. Derivation of the effective operators for the coupling configurations 

1. Z configuration 

We begin with the Z coupling configuration, which is similar to the X configuration. Both are required for the generation of 
GHZ states. We use iTint.z, F^drive.z from Eqs. ( |S18| l and ( |S21| l and the Lindblad operators in Eqs. ( jSl l] l-( [ST6l l to set up the 
non-Hermitian Hamiltonian 

+ g (pJie + bJl) . (S27) 

Here we have introduced complex detunings = uie — uiz ^ ~ *7e/2 and 6^ — ^ ~ 'iKbl‘2', where F denotes the 

particular tone of the driving field. Eurthermore we have changed into a frame rotating with the frequency of the drive ujP . Eor 
the derivation of the effective operators we will for simplicity drop the sub- and superscripts denoting the coupling configurations 
and field tones. To invert the non-Hermitian Hamiltonian we divide it into four blocks 

A = 6b^b, B =gb^Jie, C = gbJ^ D = AJee 


(528) 

(529) 
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After this separation we can formally invert the Hamiltonian using Banachiewicz’ theorem ll5^ for the blockwise inversion 
of a square matrix, 


— a b c di (S30) 

d= {D- CA-^B)~^ , a = A-i + A-^B d CA-\ b = -A-^B d, c = B'^ = -d CA-^ (S31) 

We now need to compute d to obtain any of the above elements. As we shall see, it is possible to simplify the calculation 
and to obtain closed expressions for the decay rates if we separate the involved operators by the number of atoms in state |1). 
This is done by introducing projection operators which project on the states with the same number ni (from now on, n) of 

atoms in |1). For the Z and X conhgurations discussed here, ni (or n_ for X), is conserved under the couplings by the coherent 

interactions Hint and V, but can be changed by the dissipative jump processes Lk, e.g. from ni to ni — 1 in the case of Z. Using 
these projectors we can split the non-Hermitian Hamiltonian of the excited states and its four blocks by n, 

N N 

^ ^ An + Bn + Cn + Dn (S32) 

n—0 n—0 

The inverse of these non-Hermitian Hamiltonians for each n are then found to be 

-^NH,ra = (In + bn + Cn + dn (S33) 

The effective operators of Eqs. ( |S24| i-( |S25] l are formally given by 

d/n.eff — ^ ^ b bn^Pn 

n 

-^70,a,eff — ^ ^ 'V^'TO 10)a (^| Pn 
n 

^ ^ 'V^'Tl I l)a (^| Pn 

n 

HnS = PvY,dnVPn + H.C 

n 

To obtain the effective Lindblad operators and the effective Hamiltonian it is thus sufficient to compute the blocks dn and of 
the inverse non-Hermitian Hamiltonian. Using the identities a{a^a)~^a^ = 1 and PgJi^Jee = PgJie (where Pg is the projector 
onto the ground states) we obtain 


(534) 

(535) 

(536) 

(537) 


dn = — 

A 


A5 


-1 


dee I tr o 1 d-i 




U - 9 i j 

bn — ~ Jle 

AS 


AS 


-1 


dee I „ I diddle 


(538) 

(539) 


With this, and readopting the sub- and superscripts for the conhguration and the field, and changing to a more detailed notation 
for the effective Lindblad operators we hnd 


(Fi ( 


ni 


nig 


^ ni — ^ ^ 


ni,F ^9z,n 


L. 


= 2^ -9-e 


ni ,F 


9A^ '' 
F ^^Z,ni 


L'yj^^^a,Z — 'y ^ 

ni ,F 


'JFfP'z ^ 


V / ni.f 

\ |l)a(l|F, = E pEP 

/ ni,F 


x(F) nig 

~,iF) 


(540) 

(541) 

(542) 
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To obtain this we have used the identities PgJieJeiPg = PgJiiPg, JiiPni = niPm, and JeeJeiPg = JeiPg- In the last 
expression we have also introduced the effective detunings Az,ni and the effective couplings gz,ni with 

^2 

(S43) 




nig^ 


r(^’) 


~(F) 

ffz.n, = 9 


nig 


(S44) 


As can be seen from Eqs. ( |S40| l-( [S42l l, dealing with multiple frequencies in the drive leads a priori to a sum over terms for all 
fields in the effective Lindblad operators. However, as the frequencies of these fields are well-distinguishable, we separate the 
Lindblad operators by their driving field F. Given the quadratic appearance of the Lindblad operators in the master equation, 
we can also drop the exponential phase factors. For the effective Lindblad operators for the fields F we then obtain 


N 


= E 




PF) 
■^70,a,Z 


PF) 
^71,a,Z 


ni =0 2 p 2 .il 
^ /- r,(F) 

ni—0 
^ 


ni=0 


2K^P 

Z,ni 


We also define the corresponding effective decay rates 




(F) _ ) 




(F) 

To.Z.ni = 


(F) 

Tl.Z.ni = 


4|SiP 


Z,ni 


4|A 

4|aSJ^ 


(S45) 


(S46) 


(S47) 


(S48) 


(S49) 


(S50) 


The operators in Eqs. ( |S45| )-( |ST7| ) are then the effective Lindblad operators for the Z configuration. As can be seen from the 
expressions in Eqs. ( |S43| l-( |S44l i, the effective detunings can be made very small by a suitable choice of the frequencies 


XF) 


of the fields F which can be used to engineer the rates 'yP and of the effective decay processes. The engineering 


(F) 


of the effective decay process to prepare GHZ states will be subject to Section IV 


We now turn to the effective Hamiltonian. The effective Hamiltonian is computed from Eq. ( |S24| ). Other than for the effective 
Lindblad operators, introducing a multi-tone driving field results in cross terms between different fields, here denoted by F and 
G, 


Hz = - 


N 

EE 

m^O F,G 
N 




~^{F) 

Z,nx 


nig^ 

~FF) 


= -2. A. 

ni =0 F,G 


Re 


ninPnf[^_^^^(F)_^(G)P 


4A 


(F) 

Z,ni 


./ (F) (G)x, 

3-*(‘^z -“z )‘ p 


P P 


ni + H.C. 


(S51) 


(S52) 


Here, all terms F ^ G have fast rotating exponential phase factors. Restricting the treatment to F = G where these terms 
cancel, we obtain the main contribution 


N 




1 4A^)^ , 


N 

p = V p 


(S53) 


ni =0 F \ ‘^^Z,ni / ni =0 F 

We thus find that the main effective Hamiltonian processes are AC Stark shifts with a magnitude 


-_Re 
*Z,ni — 


ni{dpf 

{F) 

Z,ni z 


4A 


As we will see further below, our choice of the field tones will make these Hamiltonian terms compensate each other. 


(S54) 
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2. X configuration 


We perform an analogous treatment for the X pumping mediated by the excited level | /) and the oscillator mode c, using the 
couplings in Eqs. ( |S19| l and ( |S22| l. We start with the non-Hermitian Hamiltonian 

= ^Pjff + + gcU.f + gcf_f (S55) 

- ( 'p'\ ( P'\ '~(P') ( P\ 

with the complex energies = u)f — ' — ijff2 and S)^ ' = ojc — ' — zKc/ 2. Carrying out the derivation in the same 

manner as above for Z, we obtain for the effective Lindblad operators 



AF) 

FjO,a,X 


L 


'yl,a,X 


N 

E 

71-—0 

N 

E 

71-—0 

N 



^yx,n_ 




2A 


(F) 

X,n. 




2A 


(F) 

X,n_ 


|0)a(-|Pn_ 

|l)a(-|Pn_, 


(556) 

(557) 

(558) 


with the effective detunings 


The effective decay rates can be written as 


A 


(F) 

X,n- 


(F) n-g^ 


- ’ - 
— ^x 


'AF) 


~{F) 

9x,n. 


= g- 


n-g 


AF) 

^X,n- 


(F) 

7o,X.n_ 


h,X,n- 


4|Cij^ 

4|a£_P 


(559) 

(560) 

(561) 

(562) 

(563) 


These operators resemble the ones for Z pumping in Eqs. ( |S45| )-( [ST7| ) if |1) is replaced by |—). The only difference is that 
spontaneous emission is still assumed to lead to the final states |0) and |1). 


C. Reduction of the dynamics to rate equations 

Using the effective operator concept, we have so far reduced the dynamics of the open system to an effective master equation 
of the ground states. This description is exact to second order perturbation theory. In addition to that, we will later achieve a 
compensation of the effective Hamiltonian, iTeff = 0, so that the remaining dynamics are purely dissipative. We can then, in 
another step, reduce the complexity of the dynamics to rate equations of the populations. This is achieved by choosing subspaces 
of interest between which the interactions present in the system do not build up coherences. These subspaces are defined by 
projection operators, e.g. Pa and Pb- Eor negligible coherences between the subspaces, we can then trace over the Liouvillian 
evolution to obtain decay rates from subspace A to subspace B 


Ta^B = ^T{PBC{PAP^PA)PB) ~ Y. P^iPBLkPAP^PALlPB) (S64) 

k 

= YY.(^f\PBLkPAP^PALlPB\^f) 
fc / 


(S65) 
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For the subspaces we will consider, the decay rate is the same for all states in the subspace. We can then calculate the decay 
rates using a single state {ipi), 


rA^B,k = \{'^f\PBLkPA\i’^)\^ (S66) 

/ 

We will use this rate equation approach to analyze the scaling of the preparation time of the protocol in Section [VT| 


D. Numerical simulations 


The evolution of the system is simulated by numerically solving the effective master equation ( |S23| l, including the effective 
operators derived above. For the evolution due to the effective master equation of the protocol we use a Trotter-like ansatz, 
simulating the evolution under the Z and X coupling in an interchanging manner, performing base transformations between the 
eigenbases of cr^ and ax in-between. To verify our findings, we compare our result to the solution of the master equation ( |S8| ) 
after truncation to one or two excitations. 


IV. ENGINEERING DISSIPATIVE MECHANISMS EOR GHZ STATE PREPARATION 


In this section we show how the effective operators derived in the previous section are engineered to prepare GHZ states ef- 
hciently. As is discussed in the main manuscript, the engineered operators are used to empty certain states and transfer the 
population to others in such a way that in the end only the target state remains. This is achieved using the driving helds in the 
coupling configurations Z and X which activate the effective decay processes derived in Sectionjl^ By choosing suitable demn- 
ings for these driving fields we can then engineer the rates of the effective decay processes. The choice of the Rabi frequencies 
of these helds is subject to Section VI where we optimize the preparation time of the protocols. 


A. Preparing |GHZ): the Z pumping 


The hrst mechanism to prepare a |GHZ) state described in the main manuscript is the Z pumping: This process transfers the 
population of all states other than \ and to |0)®^ = ^(|GHZ)-|-|GHZ_)), and thereby, to |GHZ), without affecting 

the GHZ state. We engineer this process by applying drives 12^2 in the Z conhguration with detunings ~ ~ 'P'/F9, 

where 1 < F < — 1. Here, the subscript Z+ denotes the helds with positive detunings and Z— those with negative detunings. 

If the held index F coincides with the number of atoms in 11), ni, for a certain initial state and driving held, the held is resonant 
and the effective detunings in Eqs. (|S43|l-([S44li become 


^Z±,n-i ^Z± 

^Z± 


~{F=ni) 

9z±,m =9 


^z±°z± 

ni9 



(je + Kb) 


1 Je+ Kb 

2 y/fH ’ 


(567) 

(568) 


Since we generally work in the strong coupling limit g, the above effective detunings are small compared to those for 

off-resonant driving helds with rii ^ F, 


^z±,m ± ^ 9 

(S69) 

rii — F 

(S70) 


ni 
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With the effective detunings we can then compute the effective decay rates in Eq. ( |S48| l - ( |S50| l for the processes with F = ni 
and for the off-resonant ones with F ^ ni. 


iF=ni) _ niK,l,{d-pY 

(F=ni) _ lOeiPz 

“ (7e + ’ 

T'l,Z±,ni + ^,)2 > 


(F/ni) ^ nlubjnPf 

4(ni-F)2g2- 

(F^n-i) _ Fldej^^z 
7o.z±.ni - ^(^F-myg'^' 

(F^n-i) _ FllePz 

ll,Z±,m - 4(i^_ J^^)2^2’ 


(571) 

(572) 

(573) 


From these expressions we can clearly see that the first group of rates is engineered to be strong, while the second group of rates 
is engineered to be suppressed. Using the quantities above we obtain for the effective Lindblad operators in Eqs. (|S45|l-(|ST7li. 


L 


(F) 

K,,Z± 


N-1 

E 




Pr, 


+ yJ^zPPPN, {1<F<N-1) 


r{F) 

■^70,a,Z± 


N-1 r 

E 

ni — l 
N-1 


'yO,Z±Ji |0)a(l| + ^ [ g2 


Pm + |0)a(l|PiV, {l<F<N-l) 


L\ 


■?E±«E Pm+^l[7PlN)a{l\PN, {l<F<N-l) 


(574) 

(575) 

(576) 


From these expressions we see that for 1 < ni < TV — 1 we will always find terms with ni = F to zeroth order in g~^, 
which are much larger rates than the ones with ni ^ F, which are to second order in g~^. We can therefore drop the latter for 
1 < ni < TV — 1. The terms with rii = N which in particular affect the GHZ state need, however, to be kept. Since the effective 
detunings are engineered to depend on ni, photons scattered by resonances with different ni can be distinguished. Formally this 
is justified by the exponential factors washing out interferences between terms (cf Section III i. We can therefore separate 

the terms with different ni into individual Findblad operators, each acting on a set of states with ni atoms in |1). The enhanced 
processes are then given by 


« \JPPpPm=F, {1<F<N-1) (S77) 

4S± « ^Jlp7,i\0)a{l\Pn,=F. il<F<N-l) (S78) 

« \/^FzTJN)a{MPm^F, (1 < F < TV - 1) (S79) 

The weak decay processes affecting |1)®^, and thus, |GHZ) are found to be 

L[PP\GRZ) « y4±^”^^P„,=Ar|GHZ), (1 < F < N - 1) (S80) 

«^y4±44lGHZ) + |GHZ_)), (1<F<N-1) (S81) 

4o,ti±|GHZ) « y4E4|0),(l|P„,.^PGHZ, (1 < F < TV - 1) (S82) 

« ^V^r/±4>)-(l|(|GHZ) + |GHZ_)), (1<F<N-1) (S83) 

41ti±|GHZ) « ^JJppi\l)a{l\Pn,=NPGHZ, {I < F < N - 1) (S84) 

«^y^SJ|l)-(l|(|GHZ) + |GHZ_)), (1<F<N-1) (S85) 


The result of the engineering of the effective processes for the Z configuration is thus the Z pumping by which all population 
from states with 1 < m < TV — 1 is transfeiTed to |0)®^, and thus, to |GHZ) and |GHZ_). 

Beside effective decay processes, we generally also obtain effective Hamiltonian terms as given by Eq. ( |S53| l, which we write 
as F[z = Hz+ + Hz- with 


H 


z± 


T 


N 

E ni 

4o 

ni—O ^ F^ni 




(F). 

z ) 


N 


JF) 


-Pm ^ E E ^zPnPm- 

ni—0 Fpn-i 


F — m 


(S86) 
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These terms are AC Stark shifts, the magnitude of which is obtained from Eq. (|S54[), 


*Z±.ni 


T 


4{F-ni)g 


(S87) 


(p'] ((P) (F) (F) (F^ (F') 

By our combination of red-detuned drives ^z+'> blue-detuned drives O with ~ ^z- ^z+ ~ 

( F) 

A^z- we achieve that the shifts compensate each other. 


Hz^H 


z+ 


Hz- « 0 . 


(S 88 ) 


Given that there are no Hamiltonian terms, we can turn to a description of the dynamics in terms of rate equations: 

In Section [Til C we have described the possibility to reduce the effective dynamics further to rate equations. To perform such 
a step, we identify subspaces inside which all states have the same decay rate and in-between which no significant correlations 
are built up. It is therefore important that the effective Hamiltonian is zero. We define the subspaces by the projectors Pghz = 
|GHZ)(GHZ|, Pghz_ = |GHZ_)(GHZ_|, and for 1 < rii < — 1, which contain the states with a certain number of 

atoms in |1). The decay rates can then be calculated using Eq. ( |S 66 | ). Eor resonant Z pumping {F = ni) from a subspace with 
ni atoms in | 1 ) to one with ni — 1 due to spontaneous emission to | 0 ) we obtain the rate 

^(P’=ni)^2 


AF=ni) ... I/,,. I n r(F=ni) n lM 2 , ”l70e 


(7e + 


(S89) 


The subscripts in the above and the following decay rates specify the initial subspace, the final subspace, the physical process, 
i.e. oscillator decay (k), spontaneous emission to |0) ( 70 ) or |1) ( 71 ), and the pumping process (here: Z). As opposed to the 
resonant rate (F = ni) above, the decay rates due to off-resonant fields can be written as 
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Erom these expressions follow the loss rates from |GHZ) due to Z pumping: 
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GHZ—>-GHZ_ ,7l,Z± / ^ ^ I 7i,a,Z± II ''' 26(7^ 

a^l ^ F^l 

Here we have neglected the gain of population in |GHZ) from |GHZ_). Eurthermore, we note that due to its scaling with N'^, 
loss from |GHZ) by oscillator decay should be avoided. In addition, effective oscillator decay is not useful for the Z pumping 
process. As we will see below, this is also the case for X pumping. Therefore, we will generally choose to work with Kb/c = 0 
for the GHZ protocol (a weak cooling k il may nevertheless be used to counteract heating). The overall loss rate from |GHZ) 
due to off-resonant Z pumping is then given by 
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Here, the question mark stands for any potential final state; in the scaling analysis we will typically consider the worst state 
possible. Eor the reasonable assumption of 7oe = 7ie = 7 e /2 we obtain 
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We conclude that for p 3 > 7 , k, the rates from Z pumping, ultimately leading to |0)®^, and thus to |GHZ), are much stronger 
than the loss rates from |GHZ). 

The derived rates will be used further to analyze the error and preparation time of the protocol in Section[VT| 
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B. Emptying |GHZ_): the X configuration 


|GHZ_) is emptied by X pumping, as described in the main part: Here we make use of the fact that in the eigenbasis of 
cTx = (|1)(0| + |0)(1|), consisting of |±) = (|0) ± \l))/^/2, |GHZ_) is a superposition of states with odd numbers of qubits in 
|—), whereas |GHZ) only contains states with an even number of qubits in |—), 

|GHZ_) = ((I + + H —) + I + ••• H-h) + •••) + •••) • (S97) 

72 

|GHZ) = (I + + ++) + (I + H-) + I + •••-h) + •••) + ...), (S98) 

72 


f p\ (F') 

To depump |GHZ_) without affecting |GHZ) we apply again red- and blue-detuned laser fields ^x± detunings = 

±77 g, but including only odd field indices F = 1,3, 5, {F < N), whereas for even F = 2,4, ... we use = 0. From 
Eqs. (|S59|l-(|S60ll we find for the effective detunings for resonant {F = n_) and off-resonant (F ^ n_) excitation 


- 

A±,n_ X± 


(^F) 


-2 7/ + «c), 


-{F=n_) 

9x±,n_ =9 


g(F) 

_ Ilf + Kc 


ri-g 2 y/nZ 

For the effective decay rates we find from Eqs. (|S61|)-([S63]), 
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for resonant excitation {F = ri-) and for off-resonant {F ^ nZ) excitation. With these rates and Eqs. ( |S56| l-( |S58l l we obtain 
the effective Lindblad operators 
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EEx±= E \\/^lxtnl\^)a{-\Pn.+0(^']] + ^ a{-\Pn. , (F = 1, 3, 5, ..., (F < A) ) . 

^ ^ even n_ 

(S106) 


odd n_ 


Again, we separate the effective Lindblad operators by the frequencies of the resonances, this time depending on n_. We obtain 
similar Lindblad operators as for Z pumping, with enhanced terms to zeroth order in g~^, 


= \l^xZ'JPn -. (^ = 1, 3, 5,..., (F < iV)), 
EE”xi = TESi |0)a(-|Pn_, {F = 1,3,5,..., (F < N)), 
ElExi = ^Jzix±:l\^)a{-\Pn.. [F = 1,3,5,..., (F < N)), 


(S107) 
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and suppressed terms acting on the target state 
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It can be seen that for g ^ j,k these operators make |GHZ_) (with only odd n_) decay rapidly. The losses from |GHZ) 
(with only even n_) are caused by the off-resonant drives with F 7 ^ n_. Since these terms are of second order in 5 “^, they are 
suppressed for p ^ 7, /c. 

Again, the effective Hamiltonian is compensated by the red- and blue-detuned fields, 


Hx « Hx+ + Hx- « 0. 


We can therefore describe the dynamics in terms of rates: 

With the effective operators from Eqs. (|S107|i-(fS109[) we find the decay rates from |GHZ_), 
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For the reasonable assumption 0 / 70 / = 71 / = If 1“^ the total rate is therefore given by 
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The expression in Eq. ( |S116[ ) is the total decay rate from the |GHZ_) state which is approximately given by the sum of all 
enhanced decay rates, weighted with the number of states with the same excitation. The X pumping distributes the population 
of |GHZ_) over all other states which we here denote by a question mark and later on replace by a worst-case assumption. 

The losses from |GHZ) are only caused by the off-resonant drives with F 7 ^ n_. Using Eqs. ( |S110| )-( rSl 12| l and ignoring 
negligible gain processes we obtain the loss rates 


rGHZ->-?, 70 .X± 


rGHZ->-?,7l,X± 
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Here, the binomial coefficients originate from the number of states with the same number of atoms in |—). For the total loss rate 
from |GHZ) through X pumping with both red- (X+) and blue-detuned (X—) fields we approximately find 
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(S119) 


From this expression we see that the loss terms from |GHZ) due to X pumping are of second order in g~^ and thus suppressed 
for g Beside decay out of |GHZ_) the X pumping also causes losses from states with odd n_ which have an overlap 

with states with 1 < ni < TV — 1 in the Z basis. This affects the transport from ni = N — 1 to ni = 0 by the Z pumping and 
thus to I GHZ) by imposing a loss rate 
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~ 7/ V 


N\ FjnP) 


F J 2^-1 
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We will refer to this process as “X toss” below. As we will show below this does not, however, have a significant effect on 
the scaling of the preparation time and the error, if the strength of the X process is chosen properly. We conclude that the 
combination of Z and X pumping results in the preparation of a GHZ state from any initial state. For g ^ j,k the gain rates are 
engineered to be strong, while the loss rates from |GHZ) are suppressed. The scaling of the error and the preparation time of the 
protocol are analyzed in Section[Vl| 
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V. STRONG DmVING EFFECTS 


The operator formalism used to derive effective couplings of the ground states in the system is built on the assumption of a 
perturbative drive which is much smaller than the couplings, decay, or detunings of the excited states, <C 7 , k, g, 6, A. Effects 
from saturating the excited states that could possibly slow down the preparation process are thus not included in the rates derived 
so far. On the contrary, the dependency of all decay rates on the drives c>c suggests that the optimal Rabi frequencies O of the 
drives are infinitely strong and thereby outside the perturbative regime. A proper assessment of the driving strength fl therefore 
requires the inclusion of strong driving effects which begin to play a role for O < 7 , k. 

In the following, we include such effects in the dynamics in an approximate way. Here, we discuss power broadening and 
population of the excited states and adjust the rates derived from the effective operators to account for power broadening of 
the excited states. They will later be used to analytically derive the optimal Rabi frequencies H and the scaling of the protocol 
for strong driving. In addition, in our numerical simulations we use effective operators where we have included the power 
broadening terms manually. While this treatment is not rigorous, it provides a convenient tool for rapidly determining the 
optimal parameters. These are used to simulate the effective dynamics of the system beyond the regime of weak driving and to 
match it with the simulations describing a larger Hilbert space including excitations. In addition, we take into account the effect 
of population of the excited states in our numerics. 


A. Power broadening 

We first address the effect of power broadening (or ‘line’ broadening), regarding a simple model situation: A ground state |0) is 
resonantly coupled by a field with a Rabi frequency H to an excited level |e). In total, |e) decays at a rate 7 = 70 + 7 i, where 
7 i is the decay rate into |1). We perform adiabatic elimination by setting the derivative of the density matrix to zero, p « 0. 
In the weak driving regime, where the broadening of the excited level can be neglected, this yields an excited population of 
Pee = ^'^/j'^Poo thus an effective decay rate from |0) into |1) of 7 eff = The population gain of state |1) is then 

given by 


Pll ~ 7lPee ~ 7effP00 



(S121) 


The same result is obtained when using the effective operators. Performing adiabatic elimination with a stronger drive we need 
to take into account the population of the excited level. This yields an excited population of pee = ( 7 ^ +17^). Then, as the 

coupling of the ground state | 0 ) to |e) is increased, we need to take into account the population of the coupled subspace of | 0 ) 
and |e) rather than of state |0) only. This leads to 


” (wo?L) 

In the weak driving limit H 0, this new decay rate 7 eff = 7 !!^/( 2 ( 7 ^ + 2H^)) matches the previous result. The decay rate 
for strong driving can thus obtained from the rate which was derived from the effective operators by “broadening” the natural 
line width of the excited state, 7 ^ — 7 ^ + 217^. In the limit of strong driving H 00 , the population is found in |e) with a 
probability of 1/2 such that here 7 eff approaches the constant value of 71 / 2 . Therefore, the effective decay rate from |0) to |1) 
cannot be increased beyond half of the line width of the level mediating it. 

More generally, we now seek to include power broadening in the effective operators. This is done by replacing the complex 
effective detunings, e.g. here generally denoted as Aes and peff, by “power broadened” ones. The replacements are 

made such that the rates obtained from the effective operators agree with the ones derived by adiabatic elimination. To this end, 
we make the replacements |Aeff P —|Aeff P + and |peff P l^eff P + where n is the number of atoms that can be 
excited by the drive. From our numerical simulations it turns out that the action of power broadening needs to be doubled in the 
effective operators for Z and X to achieve high accuracy between the evolution due to the effective master equation ( |S23| ) and 
the more complete master equation in Eq. ( |S 8 | l. This can be attributed to interference between the blue- and red-detuned drives. 
Indeed, considering coherent excitation of a bright state consisting of both the blue- and the red-shifted dressed state suggests 
an increase of the broadening by a factor of two and thus the replacements of the effective detunings | Aesp —| Aesp -f 2 nH^ 
and IpeffP —ISeff P + 2nH^. While these replacements are less rigorous, they are supported by our numerical simulations: 

It can be seen from Fig. 3 in the main manuscript that (1) the analytical scalings derived using the strong driving operators 
comprise upper bounds to the numerically obtained scalings and that ( 2 ) the effective operators for strong driving agree well 
with simulations of the more complete master equation. 
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Following the reasoning above, the power broadened decay rate for the Z pumping from ni to ni — 1 is found to be 


For the X depumping of |GHZ_) we have 

27 / 
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F=i,3.... ( 7 /+ Kc)^ + 2^ 1 

and for the rate of the X process when acting on other states with ni (the “X toss” rate) 
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As a consequence of the inclusion of terms for power broadening, increasing the driving strengths in the desired processes also 
increases the effect of power broadening in the desired decay rates. In the off-resonant decay rates 7 eff oc the effect of 

power broadening is, on the other hand, negligible (given that 17^ ^ g^). The detrimental rates thus still increase for a growing 
n while the desired rates saturate. This is the limiting factor for the drive which we will use to derive the possible further 
down. 


B. Population of the excited states 


An additional reduction in the fidelity for strong driving comes from the population of the excited induced by the drive. In the 
following, we investigate this effect: 

The GHZ state, through its contribution from jl)®-^, is coupled to an excited state = (|l..le) + |l..el) + ... + |e..ll))/\/N. 
Despite being off-resonance, all tones of the driving field couple to the transition from |GHZ) to the dressed states of \ijje) and 
|l)(g)iV|i), (jj-jyjjjg Strengths of and detunings Ay2 = (V^ ± '/F)g- For sach tone the excited 

population is then approximately given by 
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We also consider the excited population caused by X pumping. This requires representing the GHZ state in the X basis as in Eq. 
(|S98|l and leads to an expression 
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These expressions are included in our numerical simulations of the effective dynamics in the strong driving regime and provide 
another limitation on O. Our analytical derivation of the scaling, on the other hand, involves upper bounds to certain error 
processes and turns out to take place in a parameter regime where the population of the excited states is not significant. Therefore, 
in the analytical considerations below we leave out its effect. 


VI. SCALING ANALYSIS OF THE PREPARATION TIME AND THE ERROR OE THE PROTOCOL 


In the following, we provide an analytical study of the scaling of the GHZ scheme. We derive an expression for the preparation 
time and optimize it by the choice of the parameters. The analysis is performed both for weak driving, beginning in Section 


VIA and for strong driving, in Section VIE 


A. Optimization of the parameters for Z pumping alone 

The different schemes presented in the main manuscript have in common, that the population of a nearly exponential number of 
states is pumped to subspaces with a polynomial number of states in a number of steps linear in the size of the system. Eor GHZ 
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state preparation we have engineered strong decay processes from — 1 > ni > 1. The decay is achieved using the Z pumping 
which is only sensitive to n\ and reduces this to ni — 1, finally leading to = 0. In order to assess the performance of the 
schemes it is thus important to know the time for a concatenated process consisting of many consecutive Z pumping steps. The 
rate of each individual decay is given by (see Eq. (|S89|l): 
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a—1 k 

The average time for this decay to occur is given by the inverse decay rate, 


2ni7oe(n^^ 
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For the total time required for pumping from an ni to an rii we add the average times for the intermediate steps 
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We can also assign a total decay rate from an initial state with ni to a final state with rii, = 1/r. It is however more 

useful to use the total preparation time and to minimize it by the choice of available parameters. Here, in particular the Rabi 
frequencies ^ of individual field tones F can be chosen, as well as the tunable decay rate 7 e. 

Since the pumping occurs from a maximal ni of TV — 1 to ni — 1, the worst case preparation time for Z pumping is (using Eq. 
(|S130[)) found to be given by 


TVi = AT—1—>^ni =0 
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(S131) 


The time from ni = 1 to |GHZ) differs from that to ni = 0 by a factor of 2. This is due to the fact that only half of the 
population is pumped to |GHZ) and the other half to |GHZ_), which is continuously depumped by the X pumping discussed 
below. Therefore, on average two attempts are required so that the preparation time is doubled. 
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As anticipated already before Eq. ( S96| l, for GHZ preparation we will from now on choose the parameter values Kb = k,. = 0, 
7oe = 7ie = 7 e/ 2 , and abbreviate 7 e = 7 . In order to obtain dimensionless optimization variables, we will furthermore write 
=: Apfl (for F = 1, 2,..., — 1) with nonnegative dimensionless variables Ap. The quantity H is a dimensionful 

frequency parameter, whose size has been chosen such that (for the weak driving calculation) all have to satisfy 

i-®- that the dimensionless parameters Ap satisfy Ap £ [0,1]. In this notation, the above GHZ preparation 

time reads: 
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where we have defined the function 
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With the same abbreviations, the error rate from Z pumping alone reads, from Eq. (S96 1 : 
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Our goal for now is to find parameters {Ap} that minimize the Z-error for any given value of the GHZ preparation time 
(or, equivalently, minimize the GHZ preparation time for any fixed error value). Minimizing G({ 2 li;’}) under the constraint 
H{{Ap}) = H = const by the method of Lagrange multipliers leads to 


Af = V 
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N-l 
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> 0 such that - ^ 
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ri ^' N — F 
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(S137) 


Approximating the latter harmonic sum gives roughly rj « {log N) /H, and we would have Ap > 1 for some F (in particular, 
for F = 1) if ?7 > 1 /{N — 1), i.e. if iJ < {N — 1) log N. 

However, one can still find the optimal assignment {Ap} minimizing G({Ai?}) while obeying 0 < Ai? < 1 and FI ({Af}) < H: 
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where we defined the “ceil- 1 ” function 


W' := { 


if a: > 1 
if a: < 1 


(S139) 


and r] needs to be adjusted such that H{{Ap}) = H. Th e assignment (S138l means that Ap = 1 for F < N/{1 + l/p) and 
Ap = r]{N — F)/F for F > N/{1 + l/ty). That (S138 i is the unique optimal solution can be checked by the Karush-Kuhn- 
Tucker (KKT) conditions ll57l Section 5.3.3], using that the function G({Ai;’}) to be minimized and the constraint function 
H{{Ap}) are both strictly convex in their arguments. 

For rj G [l/(Af — 1), iV — 1] (such that the selection between the two cases in (S1391 happens at some F G [1, iV — 1]) one can 
thus compute: 
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« ploglV ^1 - + p + log (^1 - 

« p log iV + p - (p + 1) log(p + 1) 

« p log IV . 

For p = const > 0, the error in both estimates is 0(1) as > oo and thus irrelevant compared to the log N terms. In particular, 
if one wants the numerical factor H in the preparation time (S133 i to scale like log (like optical pumping), then one needs 
p > G(l), which we achieve by letting p = const as N —>■ oo. 

The error E can generally be obtained by comparing the preparation time of the desired state of the protocol and the loss rate 
out of it. Using Eqs. (|S133|l and (S135|l, the error E is found to be given by 


E := (r„j=Ar_i_>.GHz)(rGHZ->-?.7.z) = ^ G{{Ap}) H{{Ap}) . 


(S148) 


Thus, to achieve a desired error E (which may be A^-dependent, i.e. E = E{N)), one can adjust 7 appropriately. If we assume 
a relation H = a 7 in order to limit H to the weak driving regime (and where the number a may or may not depend on N, i.e. 
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FIG. 4. Four-compartment-model of the process creating |GHZ). 


a = a(N)), we can plug this back into the worst-case preparation time (8133 1 to obtain 


= %H = (8149) 

« ^ log^ N , (8150) 

V Vj VEa'^g 

where the last estimate holds for p = const as N —>■ oo and we have neglected lower-order terms in N. The prefactor is 
minimized for p = 2, leading to a minimal preparation time (for the desired error E): 
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To summarize the optimal parameter choices for the scenario considered here: 

• The Ap for 1 < F < — 1 have to be chosen as follows (cf. Eq. ( |8138[ ) with the optimal choice g = 2): 

= ^/\2{N-F)/Ey, (8152) 

i.e. Ai. = 1 for 1 < F < 2N/3, and Ap = y/2{N-E)/F for 2N/3 < F < N - 1. This leads to: 

H{{Ap}) = ^\ogN + 0{1), (8153) 

G({Aj^}) = 21ogiV + 0(1) . (8154) 


• The choice of 7 = 7 e, and consequently of 17 = 07 and 7 oe = 7 ie = 7/2, is given by: 


2V2 gVE 

3 ^/N\ogN' 
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We also set Kb = Kc = 0. This leads to: 


T„j=jv-i-).GHZ = — 7 = /— log^ N. (8156) 

2v 2 Fa‘‘g 

The appearance of in the scaling of r„j^=jv_i_>GHZ can be traced back to the fact that the noise r_ acts on each of 
the N atoms (i.e. the prefactor of N in Eq. (|896|l). 
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B. Compartment model and effective rates 


While so far we have only discussed the Z pumping, we will now make a simplified model of both processes, Z and X, that 
create the |GHZ) state. This model is more pessimistic w.r.t. the preparation time and the error treatment than the actual 
(Lindbladian) dynamics described in Section IV Nevertheless, this model will still yield a sca ling of the preparation time as 
tghz ~ log^ N)/'/E with the number of qubits N and the error E, just like Eq. (S1511, which thus shows that this is 

indeed the best achievable scaling. 

Our simplified model is shown in Fig. and consists of four compartments as we explain now. We split the 2^-dimensional 
Hilbert space into the subspaces with ni — 1,2,..., N —1, counting the number of |l)-states appearing in a computational basis 
state (as in Section|IV]l. All states with ni = iV — 1 are put into compartment 1, as they are furthest away (in pumping time) from 
the GHZ states; the states with ni = 1,..., — 2 are then put into compartment 2. The two Hilbert space dimensions belonging 

to ni = 0 and to ni = iV are split into the compartments 3 and 4 of Fig. corresponding to the |GHZ_) and the |GHZ) states, 
respectively. The desired Z pumping process creates each of those GHZ states with equal rate ^F^ out of compartment 2 where 
using Eq. ( S130| l we find 


— ('t'ni—A^—2—J-ni—o) 
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Here, we have again plugged in 7 oe = 7 / 2 , Kb = 0 and the optimal Aj from Eq. (S152 1 and followed the same computation as 
for (|S143|l, neglecting lower-order terms. Similarly, the rate from compartment 1 to compartment 2 is given by Eq. (|S128|l: 


12 


= r, 


(_F=1V-1) 

ni—N — l—¥ni—N—2,'yO,Z 




(7e + «f))^ 


7 


(S158) 


As compartment 1 is the furthest away from the desired GHZ state, we pessimistically model the X-toss (i.e. the action of the X 
process on the states other than |GHZ) and |GHZ_), see Eq. ( |S120| l) to throw any state back to the ni = N — 1 compartment, 
with a rate E^®® to be computed below. Similarly, we model the errors E^ and affecting |GHZ) such that they move the 
|GHZ)-state back to compartment 1. By the same rationale, even the “good” X process EJ^ is modelled to take the un-wanted 
state |GHZ_) back to compartment 1. 


The detrimental Z rate for the compartment model as computed with the parameters from the previous subsection, Eq. (S152i 
(see also Eqs. ( |S135 i, ( S136| l, and (S147 1 with the optimal 77 = 2) is 


Tz = rGHZ->; 


.7.Z - 


37 ( 2 ^ 


NG{{Aj}) = 


7 ( 1 ^ 3N\ogN 




8 


(S159) 


For the X rates, it was found in Section 
to Section 


VIA 


we write = A^ 


IV B 


(FJ 


Ej^ for the compartment model of Fig. |^is then given by Eq. (S116 1 : 


( p\ 

that only the Rabi oscillations 12)^-' with odd index F should be turned on. Similar 
12 with the dimensionless parameters A^S^ which we discuss below. The “good” X rate 


rt = Eghz_- 




27 / 


( 7 / + 


£ (") 


_ 212^ A 1 

' F=l,3,... 


N\ 


F{A^P)\ 


2N-1 
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Since the (normalized) probability distribution {(^) /2^ ^}f=i 3 strongly peaked around the values F ^ N/2 (similar to 

the full binomial distribution), the exact functional choice of the X coefficients A^^ (for odd F) does not really matter in the 
limit of large N, as long as it is not exponentially fine tuned. Since a similar binomial distribution occurs in the X error rates 
below as well and as the same reasoning applies, we can take all ^ to be equal as a very good approximation; 








(SI 62) 


With this, the rate Ej^ evaluates to (with exponentially good accuracy for large N): 


. 212^ N 

E+ =- At = 

^ 72 ^ 


12 = 


NA%, 
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As we will see below, with this choice the error induced by the X error rate is smaller than the one from the Z process and thus 
not a major limitation to the preparation of the entangled state. 

The X error rate is given by (S119 1 , and we simplify it again with our above parameter choices, making approximations of 
the binomial distribution and the other sum which are good in the large-7V limit; 


Tx = rGHZ-)-?,7,X 
N 


1 


7/ 

2g2 2^-1 

^ n=0,2,... 


N 

n_ 


n_ ^ F 
F=l,3,... 


n 


in 


F — 71- 


iV V 


F 


2 {F-[N/2]e,enf' 
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In the expression above, [Af/2]e^en denotes the next higher integer number of N/2. Having the limit of large N in mind, we 
evaluate the last sum as follows; 
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F odd 
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F 


F odd 

E 


F 


[^1‘Aeven) 1 <F<[nJ 2]^^,„-1 [N/2]even) [N/2],^tXl<F<N 
[7V/2]e„en , F - [iV/2], 


— + E 2 

l<F<[X/2]e„e,.-l ~ [^1‘Aeven) l<F<[tV/2]„e,.-l ~ [^1‘Aeven) 


F odd 


E 


[^/2]e 


Fodd 

E 


F-[N/21 


F [N/‘2]even) [Ar/2]„er. + l<F<Ar F [■^/2]et,en) 


=[N/2] 


[N/2\ 

( 1 - 1 

even ^2 32 


([iV/2]e„en - 1)2 

1 


1 1 

1 + 3+--- 


1 


i[N/2]even " 1) 


+ [iV/2],„en ( ^2 + 32 + • ■ • + - 1)2 


1 


. . . + ■ 


^ ^ 1 „ iV 3 1 , 3 f^ 5iV 

r:!2-— > ^ = 2- > -F = N -« 1.23iV «- 

2/Z„2 24-^ n2 4-6 4 

n=l,3,5,... n=l,2,3,4,... 

neglecting subleading terms in N. This finally gives; 


1 ' 3. {N -[N/2]even-l) 
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Finally, the X toss rate Tx'*® in Fig.|4|is given by Eq. (S120l, which with our parameter choices becomes [note that (S120i is 
half of the “good” rate (|S116[), which we have computed in Eq. (|S163|l already]; 


-ptoss 

^ X 


NA\ 

~^F~ 


(S173) 


Erom the process in Eig.|^one can see that the optimal parameters Ax have to be chosen such that the good X-process and the 
good Z-process have about the same rate; If the X pumping rate Fj is too weak, population will accumulate in |GHZ_) by the 
Z pumping (F^). On the other hand, a too strong X pumping will hinder the preparation mechanism through the X toss effect. 
We thus set the rates for the desired processes, Z pumping and X depumping to be equal. 


r+ — F+ 
^ X — ^ z- 
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This results in; 


42 - 


SJVlogN 


(S175) 


(note that, for all N >2, this choice is consistent with the requirement Ax < 1 for the weak-driving analysis). 

With these choices made, we summarize the parameters and effective rates for the four-compartment model of Fig. found 
far; 
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• For the Z-pumping, we make the choice of parameters found to be optimal in Eq. (S152 1 : 

r 1 


4(-P) _ 
— 


for F < 2N/3 , 
forF>2iV/3, 


meaning that . 


• For the X-pumping we take (see Eq. (S175 1 ): 


A 


(F=odd) 


1 


3 A^logTV 


A 


(F—even) 
X 


= 0 . 


Then one obtains for the effective rates in Fig. 


o2 2 

r+ = r+ = — 

^ ^ 7 31ogiV ’ 


J^co7np _ 


Y^toss _ 

^ X — 


Tz = 


2 = r+.31ogA^, 


1 

1 = ir+ 

7 31ogiV 2 ^ ’ 

')Vl? "iNXogN 


r 

7f22 


8 

57V 


r“ = 

^ ^2 241ogA^' 


• The total error rate (leading from compartment 4 to compartment 1) is thus 


r_ r; + rj = 7 ^ ^ 

^ ^ 9 ^ 8 \ giog^A^ 


_ 729iVlog^7V 
“ z ‘2 


9 


16 
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C. Transition matrix, stationary error, and GHZ preparation time 


The transition matrix for the 4-compartment model of Fig. is: 


/ j-\comp 

Ytoss 
^ X 

r+ 

^ X 


J^comp 

-(rr^ + r+) 

0 

0 

0 

ir+ 

2 ^ z 

ir+ 

2 ^ z 

_r+ 
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0 

V 0 

0 
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r+ 

^ z 


/—31og 
31og 
0 
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r-/r^ \ 

0 

0 
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The steady-state population poo '■= (Ai(oo), ^ 2 ( 00 ), P^{oo)^ Pi{oo)) is given as the solution (normalized to the sum of entries 
being 1) of the equation Tpoo = 0. This gives: 


Poo 


/ Pi{oo) \ 


/l/logiV\ 

P2{oo) 

_ 


A 3 ( 00 ) 

V P4{oo) / 

■ 

1 

V r+/r_ / 


The steady-state fidelity is just F = ^ 4 ( 00 ), and the error is thus 


E = I-F = l-P 4 (oo) = 1- 
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(3 + 


1 

log N 


r_/ 7^ 27iVlog^iV / 5 \ 

r+^ logA^y <?2 16 ^ giog^A^y 
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(Here the approximation in the second line was made for analytical convenience and gives a slightly pessimistic bound.) In the 
scaling with large N, this expression for E agrees with the one implied by Eq. ( |S155[ ) that was found by other means before, 
and the prefactor is similar. Thus, to achieve a desired stationary error E, we need to adjust 7 such that; 


= g^E 


277Vlog^7V 


16 


1 + 


9 log"^ N 


1 + 


1 

3iogivy_ 


- 1/2 
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Below we will use this expression instead of the results obtained in ( |S155 1 which were derived by considered only the Z pumping. 
So far we have discussed the Z pumping separately from the X pumping, deriving an individual characteristic time 
t„j= 7 V_i_>Ghz- It now remains to derive an analytical expression for the total GHZ pumping time that is obtained in the 
presence of X pumping, using the parameters ( |S188[ ) and ( |S152| l. In order to factor out the dependence of the stationary error E 
and to obtain a tractable analytical expression, we make (for the computation of the GHZ preparation time) an approximation to 
the transition matrix (|S184|) by dropping the small terms leading out of the GHZ state (these terms vanish in the limit of E —i' 0). 
That is, we set the fourth column in (|S184|i to zero: 




/-3 log TV i 1 0\ 
31ogTV -| 0 0 

0 5-10 

VO I 0 0/ 
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For the initial population vector (Pi(0), ^ 2 ( 0 ), ^ 3 ( 0 ), ^ 4 ( 0 )) = (0, 0,1,0), which corresponds to the whole population being in 
the worst state |GHZ_) of Fig. we have the following evolution: 


/Pi(f)\ /Pi(0)\ 

I Mt) _ tT+ 

Psit) - 
\PA{t)J 


P^iO) \ _ (tr+)(T+/r+) (0 
PsiO) ~ 1 

\pm/ Voy 
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Thus, on time-scales larger than any fixed to, the transition from the worst state |GHZ_) to the desired state |GHZ) happens at 
least as fast as in an exponential decay process, with an approximation F^ for the effective exponential rate computed as: 


AM = 1 u. r+= G 


r^fo 


(S191) 


Note that the fraction in the last expression will depend both on (F^fo) and on TV, since the TV-dependence of the transition 


matrix T+ in Eq. (S1891 cannot be factored out completely. To get a meaningful expression for the exponential rate, the timescale 


to should be chosen comparable to the other relevant timescales of the process. For definiteness we will thus set To := l/rj 
throughout. 

The characteristic preparation time of |GHZ) in the sense of an exponential rate is then; 
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i -(- i /7 
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In the limit of large TV, the square root inside the bracketed expression will tend to 1. And also the fractional expression involving 
P 4 (To) will tend to a constant number independent of TV, since large values of the hrst column in the matrix in (S189 1 mean that 


the transition time out of the first compartment is insignihcant compared to the other transition times, which are all independent 
of TV. This is also easily seen numerically. 


In the following table, we evaluate the factor in square brackets (which we call &(TV)) in Eq. (S194i for different values of TV: 


TV 

2 

3 

4 

5 

6 

7 

8 

10 

20 

50 

100 

500 

1000 

10^ 

10^ 

10® 

b{N) 

55 

33 

27 

25 

23 

22 

21 

20 

18 

17 

16 

15.5 

15 

14.6 

14.3 

14.1 
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The GHZ preparation time is then: 


tghz = b{N) 


VNlog^N 

a'^g\/E 


(S195) 


where b{N) tends to about 13 as iV —>■ 
and (|S188|l together with H = a'y, 7 oe 


oo (see table above). The parameter choices for this can be found in Eqs. (S176 1 , (S177 1 , 
= 7le = J, Kb = Kc = 0. 


D. Analysis of the “dynamical problem” 


Here we analyze the “dynamical problem”. This means that, for a certain exponential form of the time-evolution of the GHZ 
error E{t) (or, equivalently, of the GHZ fidelity F{t) = 1 — E{t)) in an effective model, we compute and minimize the time t 
it takes to achieve a desired target error £. The main result will be that, up to an additional factor of log(l/£’), the scaling of the 
preparation time with the particle number N and with the error S is the same as in the previous analytical approaches, see e.g. 
Eqs. ( |S156| ) and ( |S195| l (with E replaced by £). 

In the following, we write again H = 0:7 to be able to limit our treatment to the weak-driving regime by suitable choices of the 
number a. Euithermore, to keep the main derivation as general as possible, we write the rates into and out of the desired GHZ 
state as 


q2 

r = —f{N) = aVW , 

r_ = ^/i(7V) 
r 

with functions / = f{N) and h = h{N). Later, we will evaluate our results for the functions 


h{N) = 


3N log N 


1 + 


Qlog^iV 


which is motivated by Eq. (|S183 1, and 


fiN) = 


3 -I- 9 log 


which is chosen such that the ratio r_/r yields the stationary error from Eq. (S187i. 
Einally, we make the following basic ansatz for the time evolution of the error: 


m = ^ 
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(S198) 


(S199) 


(S200) 


Note that in the limit of f —> 00 the error indeed converges to r_ /E with an exponential rate given by kE, where k > 0 can be 
a dimensionless constant to adjust the effective decay rate in a model where the effective decay rate kE does not match with the 
steady state error E_/E (see Eqn. (S222 1 and below). Since we are interested in the regime of small stationary error E_/E, we 
can approximate and continue: 
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E_ 

Y 


g-tfcr ^ 


7 

9 vlW)) 


Y HN) 

7 m) 


-fexp [—jtKYf{N)] 


+ exp 


7 

9 77W)) 



mrm 

/l(7V)l/2 J 
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where we abbreviate with a constant c and a “rescaled time” r as follows: 


7 Tm) 

9 77W)' 

2/(iV)3/2 


r := tKga 
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We can treat c as a free optimization variable, since 7 is a freely adjustable parameter, and thus we can adjust c to any non¬ 
negative real number by choosing 7 appropriately (even if g and h{N) and f{N) are hxed). Furthermore, r is essentially the 
same as the physical time t, but rescaled by a fixed number (which in particular depends on N and g). 

Now the dynamical problem is as follo ws: Giv en any fixed target error 8 S (0, 1), we would like to find c > 0 such that the time 
T needed to achieve this error by Eqn. (S203 1 is minimized. Obviously, from (S203 1 , any such suitable c satisfies c S (0, y/8). 


Thus, we can explicitly solve Eqn. (S203 1 for r given c and 8\ 

-\og{£-c^) 
T = r(c) = -- 


(cG(0,V£)). 


(S206) 


To minimize this (rescaled) time r = r(c), we set its derivative equal to zero (note that a minimum exists since limc-yo t{c) = 

dT{c) 2 log(f — c^) 


t(c) = + 00 ): 


= 0 










- log(f - c^) = - 2 - 1 - 
1 


8 — 
-28 
8-c^ 
-28 
8 


= exp 


exp 


-2 + 
28 


28 

8-c^ 
28 


£-c 2 

= W{-28l^) , 


8 — (P- ^ 

= -28e-‘^ 


(S207) 

(S208) 
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where W is (a branch of) the Lambert Wfunction, which satisfies W = z {W is defined to be a solution to this equation 
for every z g C). Now we have to identify the correct branch of W{z) for our purposes: Eirst, since 8 G (0,1), the argument 


—2£/e^ in (S211 1 satisfies—2£/e^ G (—2/e^,0) C [—1/e, 0); secondly, since c G (0, v/S), the image in Eq. (S211 1 satisfies 


—281(8 — c^) G (—c», —2). Both these things together mean that the correct branch (solution) of the function W(z) in our 
problem is the branch kF_i : [—l/e,0) —>■ (— 00 ,—1], z ^ W-i{z). 

Then we can continue from (|S211[), and solve explicitly for the time-optimal c: 


= v^Jl + 


kF_i(-2£:/e2) ■ 


(S212) 


The optimal time r can be obtained by plugging this expression into ( S206| l and simplifying the expression, but there is a less 
direct and somewhat easier way: Observe from Eqs. (S206 1 and (S207 1 that 


2c C 28 [Eg. jS211| l] —W-i{ — 28/e^) [Eg. l |S212p 

£-c^ ^ 8 8-c^ ~ 8 ^ 


= VVF_i(-2£’/e2)2 + 2 VF_i(- 2 £/e 2 ) 

y8 


1 , 1 


(as 8 —)■ 0 ), 


(S213) 
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where in the last step we used the asymptotic approximation of our branch of the Lambert W function: W-i{x) = log(—cc) — 
log(— log(—a;)) -I- 0(1) as a; —)■ 0. This is justihed when we are interested in the case of very small or asymptotically vanishing 
error £ — 0 . 


Einally, we can use the above expressions to solve Eqs. ( |S204| l a nd (|S205[ ) f or the physically interesting optimal parameters 
7 = 'y{N,8) and fcHZ = t{N,8) using the values given in Eqs. (S212i and (S214i. When we use the choices for f{N) and 
h{N) given in Eqs. ( |S198| l and ( |S199| , we obtain: 
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(as N — y 00 , E —y 0), 
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Note that for small dynamical error S we have W-i 0, such that Eq. ( |S217| ) approaches the previous result in Eq. ( |S188| l. 
This is due to the fact that for small S the stationary part of the error dominates. The GHZ preparation time (defined as the time 
to reach a GHZ error of value 8) is: 


^GHZ = 


T _ 27v/3 v/]Vlog^7V 

8 k a'^g ^/8 ' 

27VNlog^N 




-T 1 + 2W_1 (-;r 


la 


1 + 


91og^ 


1/2 


1 + 


1 \ 


3/2 


SlogA^, 




a 


■gVS 


log I 


(as > 00, £ —)■ 0). 


(S218) 

(S219) 


The approximations of the Lambert W function used above, i.e. -^/l + 2/W-i —>■ 1 and + 2VE_i —)■ log(l/£), are good 

only for quite small £. Eor£ = 0.1, one should instead use the exact value VE_i(—2£/e^) = —5.27, leading to \/l + 2/lE_i = 
0.788 and + 2lT_i = 4.15, which makes that Eq. ( |S219| l is by a factor 1.8 lower than Eq. ( |S218[ ). Eor £ = 0.03 the 

corresponding value is lE_i(—2£/e^) = —6.72, leading to -^1 + 2/lE_i = 0.838 and + 2kE_i = 5.63, which makes 

that Eq. ( |S219 1 is by a factor 1.6 lower than Eq. ( |S218[ ). We can also find the relation between the stationary error E = r_/r 
(from Eq. ( |S200| ) and the dynamical error £: 

5 \ 91ogAf 


E 

~£ 


_ r_ [Eqs. | |S196|| -| |S199H 7^ 3Wl0gA'' 

^ Jv ~ ^ 

[Eq. pTTel i] ^ _ 2 


8 


1 + T—^ I ( 1 + I (S220) 

VE_i(-2£/e2) ’ 

which is independent of N. Thus, using the above values, for £ = 0.1 we get for the stationary error at the optimal parameters 
E = 0.62£, whereas for £ = 0.03 we get E = 0.70£. 

It now remains to fix the value of k (appearing in ( S218|l- ( S219|l) in an appropriate way, namely such that the model (S200l 
matches the GHZ preparation process from Sections VIB| and |VIC as well as possible. Eor this, note that the effective GHZ 
preparation rate Eqhz from Section VIC can be inferred from Eq. (S1921: 


kE = Eqhz = 




r-GHZ 7 3 log W 
Using the value of E from Eqs. ( |S196 1 and ( S199[ ), we can solve for k: 

' -l 0 g(l - P 4 (fo)) \ 


l 0 g(l - P 4 {to)) 


K = 3 




1 + 


1 

SlogiV 
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Evaluating this as in Section[VIC| we get the following table: 


N 

2 

3 

4 

5 

6 

7 

8 

10 

20 

50 

100 


0.28 

0.32 

0.34 

0.35 

0.36 

0.36 

0.37 

0.38 

0.39 

0.40 

0.41 


E. GHZ scaling analysis for strong driving 


When taking power broadening into account (see Section |V A|i, then instead of Eq. (|S128|l from the weak driving scenario, the 
favorable transition rates of the Z process are now given by Eq. (|S123|l (for ni = 1, 2,..., A^ — 1): 




2ni7oe(Hy^ 


jFnj 


(7e + K&)2 + 2 ni ( H ^^=”^^)2 7" + 2 ra ^’ 

^{F=ni) 


where we have again used the parameter values and abbreviations 7e = 7, 7oe = 7/2, = ftp as in Section 

instead of (S132i, the Z pumping time is now: 

N-l 

rVii——1—fGHZ — ‘2‘'^ni—N—l—^n-i—0 — 2 ^ .yQ 


VIA 


(S224) 
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The error rate from (S135 1 remains unchanged; 


r - 37 ,, V 

ghz^?,7,z ^ (TV - T^)2 ■ 
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As below (S136 1 , we now minimize rGHZ->?, 7 ,z (as a function of the variables {fl/}) while keeping r„j=Ar_i_>.GHZ constant. 
This leads, by the method of Lagrange multipliers, to: 


N — F 

n% = (forF = l,2,...,Af-l) 

r 


(S228) 


where A is a Lagrange multiplier (that has units of frequency^), which we will determine later. Plugging this back into (S226 1 
and ( S227| l, we get; 


N-l 


— 1 ) 27 ^ 

tn-i^ghz = -:- 1 —^ 2 ^ 


7 


F=1 


1 4(Ar-i) 27 , 

1\T -P - - - + 

N — b 7 A 
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rGHZ->?,7,Z = 


37 A 

I 652 


AT log TV. 


Thus, the stationary Z eiTor is (cf (S148 1 ): 


Ez = (rGHZ->?,7.z)(T'Ar-l->.GHz) = (A^ — 1) log + ^^A^log^A^. 


(S230) 


(S231) 


Thus, when the desired stationary Z error Ez is given, then A and 7 are determined by each other (since g and N are fixed): 


A = 


Sg^Ez — 37 ^ A^ log^ N 
6{N -1)N log N 


for 7 ^ e 


Sg'^Ez 
3N log^ 


(S232) 


Plugging this back into (S229 1 , we get; 

TAr-l->.GHZ = 4:{N — 1) 


Sg^Ez _ 2 

3N log2 N I 


(S233) 


For fixed N, g, Ez, and 7 , this is the minimal Z pumping time (i.e. minimized over all choices of pumping rates {12/}). Since 
we do not want to fix 7 a priori, we minimize the last expression over 7 , finding that the optimal choice is 


7 = 


I Sg^Ez 

9N log^ N 


The corresponding pumping strengths can now be computed from (S228 1 and (S232 1 : 
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Sg^Ez 


N - F 


9N{N -l)logN F 


(forF = l,...,TV-l)^ 


and the optimal Z pumping time for given N, g, Ez is then (compare to (S151 1 without 


TAr-l->.GHZ 


9 (AT-l)v^logA^ 

V 2 gVEz 


9 AT3/2logAT 

V 2 gVE^ 


From now on we set {N — 1) ~ AT, as this neglects only subleading terms. 
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We also take power broadening into account for the desired X pumping rate (see (|S124|i) and for the X toss rate (see |S125|l, 
which we both compute as in (|S160|l-(fS163|l and (|S173|l, and we take the detrimental X rate from (|S172|i: 
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Now we adjust the X parameters Q. and 7 / such that the X rates agree with the corresponding Z rates (cf. Section VIB 1 , i.e. 

- VI 


gVEz 


9 x372iogx and = T^hz^? = # xa/EgX - solve this for 7 / and n 


3/2 


r+ = 2r^- = r+ := l/(rAr_i^GHZl = 

exactly, one would have to solve a cubic equation. As this is quite cumbersome and uninformative, we are looking for solutions 
which satisfy the following scaling ansaetze; 7 / ~ iV“(logiV)'5 and n ~ Af‘^(logThen, one finds actually two possible 
solutions for Q and 7 /, which lead to the desired scaling of and (as N becomes large). One of the solutions is given by 
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3.51/4 /V 3 / 2 (log/V)l /2 
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(S242) 


V5 /Vi /2 ■ 

Formally, another solution exists, but this is in the extremely saturated regime where our effective operators do not apply. 

Finally, we need to find the relation between the above rate F^ := l/(TM-l-^GHZ,z} and the total GHZ preparation time 
TGHZ = l/r+ (which includes the errors and X toss) and the total stationary error E, just as in Section VIC For this, we 
consider a simplified 3-compartment model which constitutes a very good approximation in the large-A^-limit, in which the 
parameters ( |S241[ ) were computed in the first place. Here, compartment A comprises the states with ni = 1, 2,..., — 1, 

compartment B is the state |GHZ_) and compartment C the state |GHZ). Then the transition matrix is (cf. (S184 1 and Fig.[^; 
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again with F_ = F^ + F^, and again as in Section VIC we denote by T+ the transition matrix without the “bad” rates F_. 


From the stationary vector p^o satisfying Tpoo = 0 we can again compute the total stationary error E: 
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We therefore have to set Ez = E/'t) in all previous expressions. The relation between Fj and F+ = 1/tghz is similar to Eq. 

dST^ : 


r+ = 


+ -log(l - Pc{h)) 




(S245) 


where again Pc denotes the population in the GHZ state when starting from |GHZ_) state. Computing this for the above 
transition matrix and plugging in the values for F^ and E from above, we obtain; 


F+ = 0.216-F^ = 0.0339- 


g\/Ez 

Ar3/2 logJV 

Thus, in the large-W limit, the final GHZ preparation time tghz is: 
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This is achieved by the following parameter choices in terms of g, N, E: 


7 = 0.42- 
Qf = 0.42- 
n = 0.24- 


gVE 

^/N\ogN 

gVE 

gE^/^ 


(this is the 7 -rate for Z pumping), 


N-F 
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(forF = l,...,iV-l), 


iV3/2(logiV)l/2 
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Again, as in Eq. ( |S221| l, if one is interested in the dynamical error 8 instead of the static error E, one should everywhere set 
E = 0.628 for 8 — 0.1 (oi E = 0.708 for 8 = 0.03). Furthermore, the GHZ preparation time is then prolonged by an additional 
factorof log(l/f) (see Eq. (S219i). 


F. Dynamical problem for fixed preparation time T 


If the maximal preparation time T = tmax of the preparation procedure is limited, e.g. by experimental constraints or by the 
available coherence times of the underlying hardware, we can still try to adjust the remaining parameters of the preparation 
scheme in such a way, that the error 8 = E(T) = 1 — F{T) after this fixed time T is minimized. We first take the same setup 
of the dynamical problem as in Subsec tion|VID[ i.e. as in Eqs. (|S196 i-(S197 1 and (S201 1 ; 
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where we again assumed a relation = a'y between H and 7 with a fixed parameter a. It is now our task to minimize the 
preparation error E{T) for each fixed T, N, g, and for the given constants and functions a, k, f{N), and h{N). The free 
parameter is thus 7 (which determines H = 07 ). 


This optimization of the error E(T) for fixed preparation time T can be done similarly as the optimization in Subsection VID 
To find the optimal parameters 7 and fl, we demand -^E(T) = 0, which becomes: 


Tna ^= 0 . 


This is equivalent to 


s" /(W) 




whose solution is again given by the Lambert W function: 


jTmf{N) = Wo 


T^K^a'^f^{N)g^ 
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where this time we have t o take the branch Wo '■ [0, cxd) —>■ [0, c»), z 1— )• Wo(z) since both sides are nonnegative. We c an now 
solve the condition (S256 1 for the optimal f{N) pj^g back into the optimal preparation error E(T) from (S254 1 


together with the solution in terms of the Lambert W function. 

Then we obtain for the optimal preparation error E{T) after a fixed time T: 


E{T) = 


h{N) W{W + 2) 

p{N) ’ 


where W = Wo 


2h{N) 
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and the optimal rate 7 is determined by ( |S257 1 . As we have seen previously, plausible behaviours of the underlying rates F and 
r_ are given by a A^-independent f{N) ~ 1 and a linear h{N) ^ N. The expression for the optimal error E{T) at fixed time 


T and particle number N can then be evaluated numerically. Unlike in Subsection VID however, an asymptotic formula for the 


Lambert W function cannot be applied here since the optimal solution lies right in the middle of the competing polynomially 
growing (in 7 ) and exponentially decaying behaviours in ( |S254[ ). 

To obtain more analytical insight, we will thus now upper-bound the particle number N, for a given fixed preparation time and 
prescribed maximal error E. For concreteness of the exposition we will also take f{N) = 1 and h{N) = in ( |S252 i-(S254i, 
although the computation is easily generalized to more general behaviours of F and F_. If we demand E{T) < E, then from 
(|S254[) we certainly need to have 


F_ 7^ 

-^ = n\ <E 

F 

and < E . 


(S259) 
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From the last inequality, we obtain the condition 'yTtia^ > \og{l/E), i.e. the requirement I /7 < Tkc? j log(l/£) on 7 . Using 
this in the first of the two inequalities, we finally obtain; 


N < E 


^ r7~i2 2 4: 2 

< 1 nag 


V (!/£■) 
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Thus, if the maximal preparation time T is given and a certain maximal error E is to be achieved, then the maximal size of the 
GHZ state is limited to a number N of particles that is proportional to the square of the “dimensionless coherence time” gT and 
also proportional to the allowed error E (up to logarithmic factors). 


One can notice that the bound on N from (S2611 could also be simply obtained (up to logarithmic factors) by solving the main 
result (|S219|l from Subsection |VI D|for N. Thus, if we now want to take the effect of powerbroadening into account, we solve 
the main final result from Section VIE for N, namely Eq. (S247 1 , to obtain (up to logarithmic factors); 


N < — {gTf/^E^/^ . 


(S262) 


Note that this same result would be obtained from the above discussion by a choice of f{N) and h{N) satisfying 
/ii/3(At) /f{N) ~ N, i.e. for example by rates F and F_ given as F = U ^/7 and F_ = /g"^. 





























